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In this talk I will discuss an effective field theory for electron excitations in graphene, and 
show how it may be approached using the same lattice field theory methods more usually applied 
to the strong interaction between quarks and gluons. I will argue that for sufficiently strong inter¬ 
electron coupling, of the same order as that found in suspended graphene samples, there is a phase 
transition to a insulating phase, described by a quantum critical point (QCP). I will then specialise 
to the case of bilayer graphene with a biassing voltage applied between the layers, so that there is a 
non-zero density of electrons on one layer and holes on the other. This is formally very similar to 
the case of non-zero chemical potential for isopsin in QCD, and will permit numerical simulations 
probing degenerate matter (albeit in 2+ Id) in the presence of strong interactions. 

1. Relativity in Graphene 

Let’s begin with a simple tight-binding model, with a Hamiltonian describing electrons in n- 
orbitals hopping between A and B sublattices on the bipartite honeycomb lattice appropriate for a 
graphene monolayer (an excellent introduction can be found in the review by Castro Neto el al [1]): 



Figure 1: The honeycomb lattice of monolayer graphene. 
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H — ~t^ ^ b\r)a(r + Si ) + {r+ Si)b(r) 


( 1 . 1 ) 


reBi= 1 


In momentum space the Hamiltonian is rewritten 



k 


where l is the length of a CC bond. Notice how the fact that a link emerging from a site is not 
matched by another leaving at 180° results in a complex-valued kernel <&(£). If we define states 

|&±) = (v / 2) _1 [a t (^) |0),then 


{k±\H\k±) = ±{<t>{k) + <£*(*)) = ±E{k), 


(1.3) 
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so that the energy spectrum is symmetric about E — 0. 

In fact, E — 0 at precisely two independent points within the first Brillouin Zone, the so-called 
Dirac points k = K± = (0, ±47t/3 y/3l). Linearise about these points; 

^ 3 

+p) = ±v F \p y ± ip x ] + 0(p 2 ) with v F = -tl. (1.4) 

The Fermi velocity vp ~ 10 6 ms _1 in graphene. Now define modified electron operators a±(p) = 
a(K± + p) etc, and combine them into a four-spinor T 1 = (fo + ,a + ,) tr . To linear order in p 
the complex structure implicit in (1.4) results in the Dirac Hamiltonian 


H = v F ^(p)oc.p^(p), (1.5) 

P 

with 4x4 matrices obeying {a n ccj} = <5 z -y. In other words, low-energy excitations in graphene 
are massless relativistic fermions with velocity « c/300. For monolayer graphene, the number of 
relativistic flavors Nf = 2, ie 8 spinor degrees of freedom = 2 C atoms/unit cell x 2 Dirac points x 
2 spins. 

Electron-electron interactions in graphene can’t be ignored for two reasons. Firstly, Debye 
screening is suppressed due to the vanishing density of free electron states at the Dirac points. 
Secondly, the effective fine structure constant is boosted by a factor c/vp and is hence of order 
unity. The following effective field theory treats interactions in a simple way [2, 3]: 

Nf r If 

S = L / dx od 2 x(\ff a yodoW a + v F Vaf^¥a + iV\fr a Yo¥a) + ^ J dx 0 d i x{d i V) 2 . (1.6) 


The field V is an “instantaneous” Coulomb potential governed by 3d Maxwell electrodynamics 
with v p <C c, interacting with relativistic electrons moving in the plane; this is a “braneworld”. In 
the large-A/y limit the V-propagator is given by 


D(p) 


m , N f w y 1 

e2 8 (p 2 o + v 2 f \p\ 2 )G 


(1.7) 


The first term is the classical result, the second is the leading quantum correction arising from 
vacuum polarisation due to virtual electron-hole pairs. In the static limit po = 0 both contributions 
yield al/r potential. The dimensionless combination A = e 2 Nf / 16££o hvp — 1 ANf /£ parametrises 
the relative strengths of quantum vs. classical, and depends on the dielectric constant £ of the 
substrate on which the graphene sits. A is maximal for “suspended” graphene for which £ = 1. 

For sufficiently large £ 2 /£ or sufficiently small Nf it has been hypothesised that the Fock vac¬ 
uum is unstable with respect to condensation of particle-hole pairs with (xj/y/) ^ 0, spontaneously 
breaking the global U(2Nf) symmetry of (1.6) to \J(Nf)®U(Nf) and leading to a gapping of the 
electron dispersion at the Dirac points [3]. In particle physics we say a fermion mass is dynami¬ 
cally generated via a chiral condensate; in condensed matter physics the resulting phase is known 
as a Mott insulator. The transition occuring at e 2 c {Nf ) defines a QCP whose universal properties 
characterise the low-energy excitations of graphene. In the (e 2 ,Nf) plane the QCPs lie along a 
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Figure 2: Chiral order parameter (xx) resulting from (2.1) for various m at strong coupling (left), and for 
Nf = 2 (right). 


phase boundary between insulating and metallic phases, and the values of exponents such as 8 in 
the critical scaling relation 

<W)U=^ oe »»h (!-8) 

where m is an explicit symmetry-breaking mass parameter, depend on Nf. This behaviour is similar 
to that of the 3 d Thirring model, whose phase diagram has been mapped numerically over many 
years [4], and inspires a treatment based on the technically simpler action (with units vp = 1) 


5 = 



Va Yu ^ V >a + iV % To Wa + 



(1.9) 


This has the identical D(p ) in the strong-coupling limit, but the Coulombic r 1 tail is screened for 
g 2 < oo. The relation between g 2 and e 2 , A is not known a priori. 


2. Lattice Approach 

Due to the absence of small parameters at the QCP, a non-perturbative approach is needed. 
We have used a lattice model based on the staggered fermion formulation on a cubic lattice, with 
action (i = 1,... ,N) [5] 

Sum = \ Ez;-%v(l + oV x )xLfi ~ X l x rinx(l ~ iS^V x _ () )x' x _ fl +m^xX- + A L V x (2A > 

Z xpi xi x 

The fermion fields % are defined on lattice sites, and the boson fields V on the timelike links. The 
sign factors = (—l) x o+--+* M -i ensure a covariant continuum limit in the free-field limit g 2 = 0. 
The action (2.1) only has a chance of recovering the physics of (1.6) at a QCP, where for instance 
details of the underlying lattice should become irrelevant. However, (2.1) exhibits a distinct chiral 
symmetry breaking pattern U(Af)(g)U(A0 away from weak coupling, there is no guarantee 

of “taste symmetry restoration” ensuring the correct continuum symmetries with Nf = 2N. 
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Fig. 2 shows the order parameter (xx( m )) as a function of Nf in the strong coupling limit 
(left) [5], and as a function of g~ 2 for the monolayer value Nf = 2 (right) [ 6 ]. In both cases a power- 
law equation of state appropriate for a continuous phase transition has been fitted; this identifies 
the location of the QCP. The left hand plot shows a critical value Nf c = 4.8(2), proving that a QCP 
exists for finite g 2 for Nf = 2 < Nf c , and is thus potentially relevant for monolayer graphene. The 
right hand plot shows a critical coupling g~ 2 = 0.609(2) for Nf = 2 , with the critical exponent 
8 = 2.66(3), to be contrasted with 8(Nf c ) = 5.5(3). The model ( 2 . 1 ) has also been used to study 
the symmetry restoring transition at temperature T > 0, expected to be of BKT type [7]. 

These results are consistent with those obtained by simulating the braneworld gauge theory [ 8 ]. 
That work initially suggested the QCP might occur at a coupling sufficiently small that suspended 
graphene samples might lie in the insulating phase. Later simulations on a 2d honeycomb lattice 
with a more realistic inter-electron potential suggest the critical coupling is slightly too strong for a 
Mott insulator to be observed, although modification of the electron dispersion in pure suspended 
samples due to the vicinity of a QCP is not excluded [9]. 


3. Bilayer Graphene 


There is just as much interest in graphene samples formed from two honeycomb layers, in 
which the A atoms of one sheet lie directly over the B atoms of its neighbour. Naively, electron 
transport should have an effective description with Nf = 4 relativistic fermions; however in the 
absence of any further interactions the Hamiltonian (1.1) actually predicts a parabolic band struc¬ 
ture [10]. Let’s rename the original intralayer coupling along a CC bond within a layer to and 
introduce new interlayer couplings t\ between overlying AB sites and 1 3 between AB sites displaced 
horizontally a distance Z. This latter coupling results in a trigonal distortion of the parabolic band, 
breaking it into four separate cones at the Dirac point. It turns out an Nf = 4 effective field theory 
description is plausible for long-wavelength excitations with kl ^tit^/t^. 

The new ingredient in the bilayer problem is the possibility of applying a bias voltage across 
the layers, inducing a negative charge on one layer due to non-zero density of electrons, and a 
positive charge on the other due to holes. Denoting this voltage by 2fi we can introduce it with 
positive minimal coupling to fields y/ on the the upper layer and negative coupling to fields 0 on the 
lower. This is formally identical to the introduction of an isopsin chemical potential in two-flavor 
QCD. The lagrangian density in continuum notation reads [11] 




( D[V;ju] +m 

V -U 


u 1 
D[V-,-ll]-mJ 



+ -L v 2 = w^+-Lv 2 , 

2 g 2 2 g 2 


(3.1) 


where D [V ; /x] = D[V;0] ■ //}!)= l)V: Lt . The parameter m is a symmetry-breaking gap pa- 
rameter due to intralayer particle-hole pairing, and j a gap parameter due to interlayer pairing; the 
bound state formed in the latter case is known as an exciton. They act as IR regulators for the 
model; in practice numerical simulations with m = 0, j ^ 0 are perfectly feasible. The model is not 
realistic in the sense that intralayer y/-y/ interactions have the same coupling strength as interlayer 
y/ - 0 ; however, from a simulator’s point of view the key feature is the identity 


deL# = det[(D + m)^(D + m) + j 2 ] > 0 . 


(3.2) 
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Figure 3: Carrier density n c in bilayer model as a function of /I (left), and in the small p j -+ 0 limit (right). 

Even with p ^ 0 the model has no sign problem and is amenable to orthodox Monte Carlo methods. 

We have implemented a lattice version of (3.1) using N = 2 staggered fermion flavors on 32 3 
and 48 3 lattices. Since Nf = 4 is close to Nf c , it is hard to identify the location of the QCP; we 
set ( g 2 a)~ l = 0.4 1 , which we believe is just to the sub-critical side of the QCP. Simulations with 
varying p have been performed, and the following observables monitored: 


carrier density n c = 

dlnSf . - 

= <mv) (0A)0>, 

(3.3) 


dln^f . - , 


exciton condensate (T A P) = 

dj =i(w$ 

(3.4) 


d\n2f . - . 


chiral condensate (T A P) = 

dm =(YV Hh 

(3.5) 


The exciton condensate (T 11 ?) (3.4) caused by spontaneous pairing of electrons and holes between 
layers is a new feature; just like any other pairing phenomenon it is expected to induce a gap 
rendering the ground state insulating. By varying p we can check how exciton condensation varies 
with the extent of the Fermi surface (strictly a Fermi line in 2+1 d). 

Fig. 3 shows simulation results for n c \ at the left the data rise steeply with p to reach a plateau 
at pa « 0.5. There is no discernable onset p Q below which n c vanishes. The plateau marks the 
saturation region where the lattice contains one fermion (both particle and hole) per site and the 
Exclusion Principle prevents further occupation. It should be regarded as an artifact of simulating 
continuous fields on a discrete lattice. Interestingly, all other known simulable lattice models reach 
saturation at pa ~ 0(1); indeed, the result for free fermions in 2+1 d is shown as a dashed line. 
In the limit T -+ 0, p coincides with the Fermi energy Ep. Since n c is directly related to the 
Fermi momentum kp, the precocious saturation is the first hint that Ep < kp, ie. the system is 
strongly self-bound. At right the j —> 0 extrapolated data are fitted to a power law, with the result 
n c (j = 0) oc /i 3 32 ( 1 ). The density rises faster than the free-field expectation n c oc p 2 coming from 
the area of the Fermi disk. 

Fig. 4 shows the corresponding plots for (T A P). Again there is a steep rise to a peak at pa « 
0.3, followed by an even more rapid descent to zero in the saturation region. The free-held results, 

^are is needed since a priori there is no reason for spatial and temporal lattice spacings to coincide at the QCP. 
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Figure 4: Exciton condensate ( V E A T) as a function of p (left), and in the small p j —0 limit (right). 


which are much smaller and extremely sensitive to j, are shown as dashed lines. The extrapolated 
data is fitted by ( x ¥ x ¥(j = 0)) oc /i 2 - 39 ( 2 ). As we will see below, this result also suggests a weak- 
coupling approach is not applicable. 
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Figure 5: Chiral condensate ( x ¥ x ¥) vs. p (left), and the ratio ( x f /x T)/y fn c vs p for various j (right). 


The left panel of Fig. 5 shows the intralayer condensate (T A T); for small p it exceeds its free- 
field value, as might be expected from the strong interactions close to a QCP, but then rapidly falls 
to zero as p increases, so that by the time pa=03 where (T A I / ) peaks, it falls below the free-field 
value. This can be understood in terms of a competition between the two condensates; intralayer 
pairing is suppressed as £> grows, because more energy is needed to excite a particle-hole pair 
from a single layer, whereas interlayer pairing is enhanced by the larger density of states as kp 
grows. Even at ju = 0 |(W) ~ j|(W )\ peak . 

Now, (T A I / ) is simply the density of exciton pairs in the ground state. In a weakly-coupled 
BCS picture of exciton condensation, the excitons would be drawn from a shell of thickness A 

i 

around the Fermi surface, where A is the induced gap; we would then expect (T A I') oc A kp oc A, 
where the last step follows from Luttinger’s theorem relating n c to the volume enclosed within the 
Fermi surface. Thus we predict 


A(ju)°c 


TO). 

vMm) ' 


(3.6) 
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The quantity on the RHS of (3.6) is plotted as a function of p in the RH panel of Fig. 5, with 
particular attention paid to the extrapolation to both thermodynamic and j -* 0 limits. For small 
p we infer A oc ^ at a QCP no other behaviour could arise, since p is the only scale. This should 
be contrasted with other models studied using lattice simulation; in the non-renormalisable NJL 
model, A oc Auv [12], whereas in two-color QCD A = O(Aqcd) [13]. In neither case is there any 
significant dependence on p. 

4. Quasiparticle Dispersion in Bilayer Graphene 

In order to confirm the picture developed in the previous section of a system with kp > p and 
A oc p it would clearly help to have direct knowledge of these quantities from the fermion (in this 
context known as the quasiparticle) dispersion relation E(k). Here we present preliminary results 
obtained from the Euclidean propagator (*P(&)¥(&)) oc e~ E ^ t on 32 3 ; to improve momentum reso¬ 
lution partially twisted boundary conditions are used [14], so that the smallest non-zero k is 7 t/ 48<2. 
We calculated the timeslice correlators in both normal and anomalous channels according to 




Figure 6: Quasiparticle dispersion relation E(k) at pa = 0.2 for various j (left), and the fitted propagator 
amplitudes (right). 


C N (k,t) = Y,(v(0,0)v(x,t))e Ca(M) = £(V'(O>O)0(*,O) e *' X ( 4 -l) 

X X 

and fitted to the forms [15]. 

C N (k,t) =Ae~ Et +Be~ E{L '- t) - C A (k,t ) = C(e~ Et - (4.2) 

The resulting dispersion E{k) is shown for several j in the left panel of Fig. 6 . The clear dip 
marks the Fermi surface; for k < kp the correlator is dominated by forwards-propagating holes, 
and for k > kp by backwards-propagating particles, as demonstrated in the relative magnitudes of 
the amplitudes A and B in Fig. 6 (right). Note that the amplitude C for anomalous propagation, 
in which eg. an electron in one layer is absorbed by an exciton transferring its momentum to an 
electron in the other, is comparable in magnitude to those for normal propagation for a range of 
momenta with k ^ kp. 
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From the plots we estimate kpa « 0.4 > jia = 0.2, supporting the earlier contention that kp > 
Ep due to the self-bound nature of the system near the QCP. As a sanity check we note that for free 
lattice fermions the carrier density n{ reelatt (jia = 0A)a 2 « 0.06, to be compared with the directly 
measured value 0.09 (see Fig. 3). 
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Figure 7: Gap A (j) estimated from the minimum of E(k). 

Finally in Fig. 7 we plot the gap A estimated from the minima of E (k) in Fig. 6, with data taken 
at two different values of /I. We can then extrapolate to j — 0 and estimate A a(\ia = 0.1) ^0.1, 
Aa(\i = 0.2) « 0.3. This is roughly consistent with our earlier prediction that A oc /j, and supports 
a physical picture of a gapped Fermi surface with A//ll ~ 0(1). This should be compared with 
studies using a diagrammatic approach predicting A/jLt ~ 10 -7 [16]. 

5. Summary 

We’ve discussed how the effective description of graphene in terms of relativistic field the¬ 
ory can be studied non-perturbatively using orthodox lattice simulation techniques familiar from 
particle physics. The simulations have confirmed at least the theoretical possibility of a QCP in 
monolayer graphene which may have relevance for the correct description of charge transport in 
suspended samples. In the case of bilayer graphene the model presented here, though arguably not 
very realistic due to the equality of intra- and inter-layer interaction stengths, is at least an inter¬ 
esting new member of the rather exclusive stable of models which permit study by Monte Carlo 
techniques in the presence of a chemical potential. As we have seen, its behaviour is rather differ¬ 
ent from other models in this class (NJL, QC 2 D) because residual interactions at the Fermi surface 
are strong in the vicinity of a QCP. 
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